function [QW_Mat,Coef_Mat]=DistApp_QW_Moment2QWCoef(QQ_Unit,CentMomOrderTable,Mean_Mat,CentMoment_Mat,Coef_Mat_Initial)

%% Preliminary
% # of discrete dimension
N_Discrete  =   size(CentMoment_Mat,2);
% # of continuous dimension
N_Continuous=   size(QQ_Unit,2);
N_QQ        =   size(QQ_Unit,1);
% # of Moments
N_CentMom   =   size(CentMomOrderTable,1);
N_Mom       =   N_CentMom+N_Continuous;
if size(CentMomOrderTable,2)~=N_Continuous
    error('Order Table does not match the continuous dimensions');
end

%% Compute and Collect the QW
if nargin<=4
    Coef_Mat_Initial    =   zeros(1+N_Mom,N_Discrete);
end
QW_Mat      =   zeros(N_QQ,N_Discrete);
Coef_Mat    =   zeros(1+N_Mom,N_Discrete);
for ii=1:N_Discrete
    [QW_Mat(:,ii),Coef_Mat(:,ii)] ...
                =   Moment2Coef_Unit(QQ_Unit,CentMomOrderTable,...
                                     Mean_Mat(:,ii),CentMoment_Mat(:,ii),...
                                     Coef_Mat_Initial(:,ii));
end

return

function [QW,Coef]=Moment2Coef_Unit(QQ_Unit,CentMomOrderTable,Mean,CentMoment,Coef_Initial)

QQ_Dim          =   size(QQ_Unit,2);
QQ_Num          =   size(QQ_Unit,1);
Mean_Target     =   Mean;
CentMoment_Target ...
                =   [zeros(QQ_Dim,1);CentMoment];
CentMoment_Num  =   size(CentMomOrderTable,1);

MaxOrder        =   max(max(CentMomOrderTable));
UnitCentMoment  =   zeros(QQ_Num,MaxOrder+1,QQ_Dim);
CentQQ          =   bsxfun(@minus,QQ_Unit,Mean_Target');
for i=1:QQ_Dim
    UnitCentMoment(:,:,i)   =   bsxfun(@power,CentQQ(:,i),(0:MaxOrder));
end

CentMomentTable =   ones(QQ_Num,CentMoment_Num);
for i=1:CentMoment_Num
    for j=1:QQ_Dim
        Order_Temp      =   CentMomOrderTable(i,j)+1;
        CentMomentTable(:,i) ...
                        =   CentMomentTable(:,i).* ...
                            UnitCentMoment(:,Order_Temp,j);
    end
end
PolyTerm        =   bsxfun(@minus,[CentQQ,CentMomentTable],CentMoment_Target');

% TempFun         =   @(PolyCoef) PolyTerm'*exp(PolyTerm*PolyCoef);
% TempFun_Der     =   @(PolyCoef) PolyTerm'* ...
%                                 bsxfun(@times,exp(PolyTerm*PolyCoef),PolyTerm);
TempQW_Fun      =   @(PolyCoef) exp(PolyTerm*PolyCoef);


if nargin>4
    PolyCoef    =   Coef_Initial(2:end);
else
    PolyCoef    =   zeros(CentMoment_Num+QQ_Dim,1);
end

% PolyCoef    =   fsolve(TempFun,PolyCoef); 

QW          =   TempQW_Fun(PolyCoef);
QW          =   QW/sum(QW);
Temp        =   PolyTerm'*QW;
Diff        =   max(abs(Temp));
ItNum           =   0;
MaxItNum        =   40;
while Diff>1e-6 && ItNum<=MaxItNum
    Temp_Der=   PolyTerm'*bsxfun(@times,QW,PolyTerm);
    PolyCoef=   PolyCoef-Temp_Der\Temp;
    QW      =   TempQW_Fun(PolyCoef);
    QW      =   QW/sum(QW);
    Temp    =   PolyTerm'*QW;
    Diff    =   max(abs(Temp));
    ItNum   =   ItNum+1;
end

M0              =   1/sum(exp(PolyTerm*PolyCoef));

Coef            =   [M0;PolyCoef];
QW              =   M0*exp(PolyTerm*PolyCoef);

return